#include "FiniteDiff/LevelOp.h"
#include "Core/TensorSlice.h"

template <int Dim>
template <class T>
void LevelOp<Dim>::LAPC_4TH(Real alpha, const Tensor<Real, Dim> &phi, const T &rhs, Tensor<Real, Dim> &LofPhi) const
{
  auto dx = rd.spacing();
  Box<Dim> valid = rd;
  assert(LofPhi.box().contain(valid));
  ddfor(valid, [&](Vec<int, Dim> idx) {
    Real r = 0.0;
    for(int d = 0; d < Dim; ++d) {
      auto u = Vec<int, Dim>::unit(d);
      r += (-phi(idx-u*2) + 16*phi(idx-u*1) - 30*phi(idx) + 16*phi(idx+u*1) - phi(idx+u*2)) / (12*dx[d]*dx[d]);
    }
    LofPhi(idx) = alpha * r + rhs(idx);
  });
}

template <int Dim>
void LevelOp<Dim>::LAPC_4TH(const Tensor<Real, Dim> &phi, Tensor<Real, Dim> &LofPhi) const
{
  LAPC_4TH(1.0, phi, [](Vec<int,Dim>) { return 0.0; }, LofPhi);
}

template <int Dim>
void LevelOp<Dim>::LAPC_4TH(const Tensor<Real, Dim> &phi, const Tensor<Real, Dim> &rhs, Tensor<Real, Dim> &rsd) const
{
  LAPC_4TH(-1.0, phi, rhs, rsd);
}

template <int Dim>
void LevelOp<Dim>::RELAX_4TH(const Tensor<Real, Dim> &phi,
                             const Tensor<Real, Dim> &rhs,
                             Tensor<Real, Dim> &JofPhi,
                             Real w) const
{
  auto dx = rd.spacing();
  Box<Dim> valid = rd;
  assert(JofPhi.box().contain(valid));
  ddfor(valid, [&](Vec<int, Dim> idx) {
    Real r[2] = {0.0};
    for(int d = 0; d < Dim; ++d) {
      auto u = Vec<int, Dim>::unit(d);
      r[0] += (-phi(idx-u*2) + 16*phi(idx-u*1) + 16*phi(idx+u*1) - phi(idx+u*2)) / (12*dx[d]*dx[d]);
      r[1] += 1.0 / (dx[d]*dx[d]);
    }
    JofPhi(idx) = w * (r[0] - rhs(idx)) / (r[1] * 30.0/12) + (1-w) * phi(idx);
  });
}

template <int Dim>
Real LevelOp<Dim>::NRM(const Tensor<Real, Dim> &lhs, int nt) const
{
  Box<Dim> bx = rd;
  if(nt == 0)
    return norm(lhs.slice(bx), 0);
  else if(nt == 1)
    return norm(lhs.slice(bx), 1) * prod(rd.spacing());
  else if(nt == 2)
    return norm(lhs.slice(bx), 2) * sqrt(prod(rd.spacing()));
  return 0.0;
}

//================================================================================
template class LevelOp<SpaceDim>;
